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Abstract. The Galactic bar causes a characteristic splitting of the disc phase space into regular and chaotic 
orbit regions which is shown to play an important role in shaping the stellar velocity distribution in the Solar 
neighbourhood. A detailed orbital analysis within an analytical 2D rotating barred potential reveals that this 
splitting is mainly dictated by the value of the Hamiltonian H and the bar induced resonances. In the u — v 
velocity plane at fixed space position, the contours of constant H are circles centred on the local solid rotation 
velocity of the bar frame and of radius increasing with H. For reasonable bar strengths, the contour H — H12 
corresponding to the effective potential at the Lagrangian points Z/1/2 marks the average transition from regular 
to chaotic motion, with the majority of orbits being chaotic at H > H12 . On top of this, the resonances generate 
an alternation of regular and chaotic orbit arcs opened towards lower angular momentum and asymmetric in u 
for space positions away from the principal axes of the bar. Test particle simulations of exponential discs in the 
same potential and a more realistic high-resolution 3D A'^-body simulation reveal how the decoupled evolution 
of the distribution function in the two kind of regions and the process of chaotic mixing lead to overdensities in 
the H ^ H12 chaotic part of the disc velocity distributions outside corotation. In particular, for realistic space 
positions of the Sun near or slightly beyond the outer Lindblad resonance and if u is defined positive towards the 
anti-centre, the eccentric quasi-periodic orbits trapped around the stable a;i(l) orbits - i.e. the bar-aligned closed 
orbits which asymptotically become circular at larger distances - produce a broad u ^ regular arc in velocity 
space extending within the H > H12 zone, whereas the corresponding <; region appears as an overdensity of 
chaotic orbits forced to avoid that arc. This chaotic overdensity provides an orig inal i nterpretation, distinct from 



the anti-bar elongated quasi-periodic orbit interpretation proposed by Dehnen (200C), for the prominent stream 
of high asymmetric drift and predominantly outward moving stars clearly emerging from the Hipparcos data. 
However, the most appropriate interpretation for this stream remains uncertain. The effects of spiral arms and of 
molecular clouds are also briefly discussed within this context. 
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numerical 



1. Introduction 

The kinematics of disc stars in the Solar neighbourhood 
displays several long known properties, such as the in- 
crease of velocity dispersion with age, the tendency of 
young stars to appear in moving groups or streams, and 
the classical vertex deviation affecting stars with asym- 
metric drift down to ~ 25 kms~^ relative to the Sun and 
mainly owing to the Hyades and Sirius streams. Disc heat- 
ing is traditionally attributed to the diffusion of stars by 
transient spiral arms or by massive compact objects like 
molecular clouds, the streams to dissolving ensembles of 
stars born at the same place, and the vertex deviation to 
local gravitational perturbations like spiral arms or local 
departures from a steady state. 



Beside these properties, the local disc velocity distri- 
bution also betrays a broad stream of low angular mo- 
mentum and mainly outward moving stars with a mean 
heliocentric asymmetric drift s « 45 kms~^, i.e. typical of 



the thick disk (Gilmore et al. 1989), which hereafter will 
be referred to as the "Hercules" stream, according to the 
comoving Eggen group C, Herculis (Skuljan et al. 1999). 



The mean outward motion of stars with high asymmet- 
ric drift, also known as the "m- anomaly" and seen up to 
over s = 100 kms^^ in metal rich samples (Raboud et al. 



1998) and in Mira variables with period between 145 and 
200 days (Feast & Whitelock 2000), is alr eady apparent 



in ea rly ste llar kinematical samples (Eggen 1966 ; WooUey 
et al. 1970 ) and was recognised long ago by Mayor ( 1972 ), 
but the clearest evidence for the Hercules stream comes 
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Fig. 1. Heliocentric velocity distribution in the u — v plane of all the Hipparcos single stars with a{n)/Ti < 0.1, d < 100 pc and 
radial velocities in the Hipparcos Input Catalogue (left) and of the sub- samples with B — V > 0.6 (middle) and B ~ V < 0.4 
(right). For the sake of comparison, the contours are as in Dehnen ( 1998| ), containing 2, 6, 12, 21, 33, 50, 68, 80, 90, 95, 99 and 
99.9 percent of all stars. The diagram for the full sample is exactly the same as in Fux (2000), except for a different labelling 
of the contours. 



from the Hipparcos proper motions combined with (Fig. |l]) 
or without (Dehnen 199S) available radial velocities. This 
stream is very likely to have a dynamical orig in be cause 
its stars are older than ~ 2 Gyr (Caloi et al. 1999) an d 
present a wide range of metallicities (Raboud et al. 1998|) . 

The existence of the Hercules stream is most proba- 
bly related to the influence of the Galactic bar. It is now 
indeed widely accepted that the Milky Way is a barred 
galaxy, as are the majority of disc galaxies. Evidence for 
the bar comes from longitudinal asymmetry in the bulge 
surface photometry (e.g. Blitz & Spergel 1991; Binney et 
al. 1997), star counts (e.g. Nakada et al. 1991; Nikolaev & 
Weinberg 1997; Stanek et al. 1997), interpretation of the 
obser ved g as kinematics in the c entra l few kpc (Binney 
et al. |l99l| ,' Englm aier & Gerhard [l999| ; Fux |l99S| ; Weiner 
& Sellwood |1999|) , large microlensing opti cal d epths to- 
wards the Galactic bulge (Paczynski et al. 1994; Kuijken 
1997| ; Gyuk |l999[ Alcock et al. ^000|) a nd possibly inner 



stella r kinematics (Sevenster et al. |1999| ; see also Gerhard 



199£ for a recent review) . Although still not very well con- 



strained, the most quoted values for the main bar param- 
eters are an in-plane inclination angle with respect to the 
Galactic centre direction (p « 15° — 45°, with the near side 
of the bar in the first Galactic quadrant, and a corotation 
radius Rcr « 3.5 — 5 kpc. 

Barred A^-body models of the Milky Way produce a 
mean outward motion of disc particles at realisti c po- 
sitions of the Sun relative to the bar (Fux et al. 1995 ; 
Raboud et al. 1998 ), but the precise bar induced dynami- 
cal process leading to the observed kinematical properties 
of the Hercules stream is still a matter of debate. Dehnen 



( |l999b| , pOOq - hereafter D2000) relates this stream and 
the main mode of high angular momentum stars in the 
observed velocity distribution to the coexistence near the 
outer Lindblad resonance (OLR) of two distinct types of 
periodic orbits replacing the circular orbit close to the 
OLR in a rotating bar red po tential, i.e. the same idea in- 
troduced by Kalnajs (1991) to explain the Hyades and 
Sirius streams. Linear theory indeed predicts that the 
orientation of orbits closing in the bar rotating frame 
changes across the main resonances associated with the 
bar (Binney & Tremaine 1987). In particular, periodic or- 
bits outside and inside the OLR radius are respectively 
elongated along the major and minor axis of the bar, and 
both types of orbits, as well as the quasi-periodic orbits 
trapped around these orbits, can overlap in space near the 
OLR. According to D2000, the Hercules stream and the 
main velocity mode, respectively "OLR" and "LSR" mode 
in his terminology, result from the anti-bar and bar elon- 
gated orbits respectively, and the valley between the two 
modes from off-scattered stars on unstable OLR orbits. 



Raboud et al. (1998), on the other hand, suggest that the 
Hercules stream involves stars merely on chaotic orbits 
and susceptible to cross the corotation radius and wan- 
der throughout the Galaxy, but do not explicitly justify 
why such stars should move outwards on the average in 
the Solar neighbourhood. One motivation for this inter- 
pretation is that of order 10% of the particles in A^-body 
models of barred g alaxie s indeed follow such orbits (e.g. 
Pfenniger & FriedH |l99lD . 

This paper investigates how the barred potential of 
the Milky Way divides the phase space of the stellar disc 
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into regions of regular and chaotic motion and how this 
segregation may explain some properties of the observed 
local stellar kinematics and in particular help to clarify 
the real nature of the Hercules stream. The investigation 
is first performed in details using the same analytical two- 
dimensional rotating barred potential as in D2000 and 
then complemented with the results from a more realistic 
high-resolution three-dimensional A^-body simulation. 

The structure of the paper is as follow: Section || briefly 
presents the observed stellar velocity distribution in the 
Solar neighbourhood and some further informations about 
the Hercules stream. Section ^ recalls a dynamical classifi- 
cation of orbits in rotating barred potentials based on the 
Jacobi integral and determines the location in local veloc- 
ity space of the class of orbits that may cross the coro- 
tation radius. Section || describes the analytical barred 
potential adopted in the 2D study and Sect. || the main 
periodic orbits supported by this potential outside corota- 
tion. Section ^ derives the associated regular and chaotic 
regions in velocity space as a function of space position 
relative to the bar. Section]^ presents the velocity distribu- 
tions at the same space positions resulting from test par- 
ticle simulations and examines the role of chaos in shap- 
ing these distributions. Section || shows how the derived 
velocity distributions depend on the initial conditions of 
the simulations and Sect. ^ how the particles initially on 
OLR orbits eventually contribute to these distributions. 
Section gives the results inferred from the 3D iV-body 
simulations. Section |ll| makes a quantitative comparison 
of the model velocity distributions with the observed one 
and discusses the most likely origin of the Hercules stream. 
Finally, Sect. |l^ sums up. 

2. Observed local velocity distribution 

There are several attempts to recover the velocity dis- 
tribution of stars in the Solar neighbourhood from the 
Hipparcos data published in the recent literature (e.g. 
Dehnen 1998 and Sk uljan et al. f 999 for al l stellar types; 
Chereul et al. 1998| and Asiain et al. 199£ for early- type 



stars). The main features of the u ~ v distribution are il- 
lustrated in Figure |l| and the mean velocities of the high- 
lighted streams are listed in Table |l|. Throughout this 
paper, v and u respectively stand for the azimuthal and 
radial velocity components, with positive values towards 
galactic rotation and towards the Galactic anti-centre, and 
w for the vertical velocity component. 

The main sample selected for this figure is built from 
the 3481 single stars of the Hipparcos Catalogue with rela- 
tive errors on parallaxes less than 10%, distances less than 
100 pc, and given radial velocities in the Hipparcos Input 
Catalogue. Here, an entry of the Hipparcos Catalogue is 
considered as a single star if the CCDM identifier and 
Multiple System Annex flag (fields H55 and H59 respec- 
tively) are void, the number of components (field H58) 
is 1 and the solution quality (field H61) is different from 
'S'. Two disjoint sub-samples are isolated from this main 
sample, the first one restricted to the 1707 stars with 



Table 1. Mean heliocentric velocities of some stellar streams 
in the Solar neighbourhood. The velocity components are es- 
timated from the left frame in Fig. |l|, except those for the 
Arcturus stream which refer to figure 3 of Dehnen ). 



Stream 


V [kms 


u [kms ^] 


Coma Berenices 


-4 


3 


Sirius 


3 


-9 


Hyades 


-18 


42 


Pleiades 


-19 


13 


Hercules 


-45 


35 


Arcturus 


-110 


16 



B — V > 0.6, representing stars which are older on the av- 
erage than the stars in the full sample, and the second one 
to the 855 stars with B — V < 0.4, representing essentially 
main sequence stars which are younger than 2 Gyr. The 
diagrams are derived using the adaptative kernel method 
described in Skuljan et al. (1999), with an average smooth- 
ing length h = 16 kms^^ for the B — V > 0.6 sub-sample, 
and h = 10 kms"^ for the other samples. 

The reader should be warned that stellar samples built 
this way are kinematically biased in the sense that ra- 
dial velocities are predominantly known for high proper- 
motion stars (Binney et al. 1997; Skuljan et al. 1999). 



Moreover, the completeness of the Hipparcos Catalogue 
depends on Galactic latitude, so that the effects of such a 
bias are even further complicated by the anisotropic local 
velocity distribution. Nevertheless, the resulting velocity 
distributions closely resemble the a sserted unbiased distri- 
butions derived by Dehnen (1998), suggesting that kine- 
matical biases do not severely affect our diagrams. 

Figure |l] nicely confirms that the Hercules stream in- 
volves merely old disc stars. According to D2000, roughly 
15% of the Hipparcos stars with B — V > 0.6 belong 
to this stream, but this is likely an underestimate of the 
corresponding fraction among local old disc stars because 
such a colour range is still contaminated by young stars 
which contribute negligibly to the stream and because the 
Hipparcos catalogue is biased towards young stars. The 
average luminosity of stars in this catalogue indeed in- 
creases with distance and the catalogue essentially covers 
the vertical region of the Galactic plane where the fraction 
of young stars is largest. 

3. Effective potential and Jacobi integral 

In a rigid potential ^(x) rotating at a constant frequency 
f2p about the z-axis, the Hamiltonian of a test particle 
expressed in the rotating frame writes: 



H{x, x) 



1 -2 
-X 

2 



*off(a;). 



(1) 



where ^cs{x) — ^{x) — ^Qp{x^+y^) is the effective poten- 
tial. If ^(a;) is non-axisymmetric and f7p ^ 0, the energy 
E and the z-component of the angular momentum are 
not conserved individually, and the only known classical 
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Fig. 2. a) Efltective potential in the x — y plane of a barred 
disc model corresponding to Eq. (j^ with F = 0.10. The bar 
is along the j/-axis and the spacing between the contours in- 
creases by a factor 1.2 towards lower ^cff. The crosses indicate 
the Lagrangian points L1/2 (on the j/-axis) and I/4/5 (on the 
a;-axis). b) Effective potential along the y-axis (thick line) and 
the three main classes of orbits related to the conservation of 
the Jacobi integral. The shaded area below the curve is forbid- 
den for orbits with H < H\2. 



integral of motion generally is the value of the Hamiltonian 
H = E — flpLz, known as the Jacobi integral. Since 
must be positive, this integral restricts the motion of a 
particle to the space region where $off (^c) < H. 

In a rotating barred potential, the contours of effec- 
tive potential in the plane of symmetry z = look like 
a volcano with a sinusoidal crest, the extrema of which 
defining the locations of the Lagrangian points ii/2 and 
L4/5, corresponding respectively to the saddle points and 
maxima of $cff on the major and minor axis of the bar 
(Fig. ^). Two critical values of the Hamiltonian are as- 
sociated with stars corotating at these points, namely 
H12 = ^cs{Li/2) and _ff45 = $off (-^4/5)- The first of them 
can be used to classify stellar o rbits into three dynamical 
categ ories (Sparke & Sellwood 1987 ; Pfenniger & Friedli 
1991 ): the bar orbits and disc orbits with H < H12, which 



cannot cross the H12 contour and are therefore confined 
inside and outside corotation respectively, and the hot or- 
bits with H > H12, which are susceptible to cross the 
corotation barrier and explore all space except a small 
region around ii H < H45 (Fig. ||b). Stars with 

H12 < H < H45 cannot cross the corotation radius at all 
azimuth and may therefore more likely be locked during 
several orbital periods on either side of corotation. 

In the Solar neighbourhood, located confidently be- 
yond corotation, only stars from the disc and hot popula- 
tions are observed. Since these stars share about the same 
$off if not too far from the Galactic plane, their iJ-values 
depend mainly on the velocities and thus one expects that 
the two populations occupy different regions in local ve- 
locity space. If V, u and w are measured with respect to 
the Galactic centre, Eq. (n^) transforms into: 



{v - RoVLpf + + ^ 



(2) 



where Ro is the galactocentric distance of the Sun and 
the local effective potential. Thus the contours of con- 
stant Hamiltonian in velocity space are spheres centred 
on {v,u,w) 
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Fig. 3. Contours of constant Hamiltonian in the u — v plane 
of a realistic 2D barred model (Eq. (j^ with Rcr/Ro = 4.5/8, 
Vo — 200 kms^^ and F — 0.20). The velocities are relative 
to an inertial frame. The inner and outer solid circles give the 
H12 and contours respectively, and the dotted circle is the 
axisymmetric limit of these contours. 



increasing with H. Stars on disc and hot orbits are re- 
spectively those inside and outside the H12 sphere. If the 
vertical dimension is neglected, these spheres become cir- 
cles with the same properties. In the axisymmetric limit 
and for a flat rotation curve of circular velocity Vo , the ra- 
dius CcR of the H12 — -^45 contour and the low azimuthal 
velocity Aw relative to Vo at which this contour crosses 
the u — axis are then given by: 



Av 



2[i/l2 -$cff(i?o)] 
2 

In ' ' 



Ro 



Ro 



Ro 
Rci 



CCR, 



- 1 



(3) 
(4) 



where i?cR — Vo/flp is the corotation radius (see Fig. ||). 
For Rck/Ro — 4.5/8 and Vo — 200 kms~\ one gets 
—Av = 0.227wo ~ 45 kms^^, which coincides with 
the mean heliocentric asymmetric drift of the Hercules 
stream^. Note however that this simple approximation is 
not truly a lower limit to the asymmetric drift of stars 
on hot orbits for several reasons: a non-zero w velocity 
component defines two circles on the H12 sphere with a 
reduced projected radius on the u — v plane (small effect, 
of order 1 kms~^ for \w\ = 20 kms~^), the presence of a 
bar lowers the effective potential at Li^2 (larger effect, of 
order 5—10 kms""'^), and finally |Au| is smaller if m 7^ 0. 
Hence most stars in the Hercules stream are likely to fall 
outside the H12 sphere and therefore may belong to the 
hot population. 

From Eqs. (^ and (^, it also follows that the radius 
CcR and the velocity separation |Au| increase for larger 
galactocentric distances relative to Ren- In particular, 
whatever the strength of the bar, one can always increase 
the fraction of the Hercules stream falling in the hot orbit 
region by reducing the value of Ro/Rcn- 



(i?of^p,0,0) and of radius ^/2(ir- 



^ It is implicitly assumed here that the azimuthal velocity of 
the Sun is close to the circular velocity of the axisymmetric 
part of the Galactic potential. This is probably correct within 
5 — 10 kms~^, as will be argued in Sect. |l^. 
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4. Working potential 

The analytical 2D barred potential adopted for the orbital 
structure analysis and the test particle simulations is the 
same as in D2000: 



2/ 1 



3/1 



5/1 



6/1 



$(i?,</)) =$o(i?)+*b(i?,< 



with 



vl In i?, 
1 



$o(i?) 



2 -Mi 



R<a 
R> a. 



(5) 

(6) 
(7) 



This represents the sum of an axisymmetric potential <I>o 
with constant circular velocity Vo and a barred poten- 
tial <i>b falling off as a quadrupole at R > a. The inner 
and outer parts of the latter component are described by 
two distinct functions which connect together at i? = a 
such as to ensure continuous potential and forces. The 
bar major axis is taken to coincide with the y-axis, con- 
trary to the convention in D2000. The parameter F is the 
bar strength, defined as the maximum azimuthal force on 
the circle of radius a divided by the radial force of the 
axisymmetric part of the potential at the same radius 
(in absolute value). It is related to Dehnen's parameter 
a by F = 8.89a. The potential is rotating at a con- 
stant pattern speed fip such as to place corotation at 
RcR = 1.25a, in agreement with numerical simulations 
and analyses of observations in early-type barred galax- 
ies if a is associ ated with the bar semi-major axis (e.g. 
Elmegreen 1996 ). Unlike D2000, only flat rotation curve 
models will be examined. In this case, the OLR and coro- 
tation radii are related via i?oLR = (1 + 1/\/2)-Rcri and 
a value of Ro/ Rqlr = 1.1 corresponds to a corotation ra- 
dius i?cH. = 4.26 kpc if i?o = 8 kpc. Some considerations 
in the case of a non-constant rotation curve can be found 
in Sect. |ll|. 

The next sections present a study of the periodic or- 
bits outside corotation in the adopted rotating barred po- 
tential, identify the regular and chaotic regions in phase 
space associated with this potential, and discuss how stars 
may populate the available orbits. All orbits in these sec- 
tions are integrated in double precision using an 8 order 
Runge-Kutta-Fehlberg algorithm (Fehlberg 1968). Two 
values of the bar strength will be considered, F — 0.10 
and F — 0.20. The larger value corresponds to a rather 
strong bar (see Sect. ^ for a quantification with respect to 
real galaxies), but has the advantage to clearly point out 
the effect of chaos in the test particle simulations. Some 
of the key results for the intermediate case F = 0.15 are 



presented in Fux (2001a). For comparison, Dehnen's sim- 
ulations were done in the range F — 0.062 — 0.116. 

5. Periodic orbits 

Much of the orbital structure in a system can be assessed 
from the study of its periodic orbits, which sometimes 
are considered as the skeleton of the system. While pe- 
riodic orbits have been widely investigated in 2D and 3D 




Fig. 4. Distinct looped periodic orbits in the axisymmetric 
potential $o of Eq. which close after one rotation in the 
rotating frame and are symmetric with respect to the x- and/or 
y-axes. The orbits all have {H — //olr)/uo ~ 0.244 and are 
not drawn at the same relative scale. For cross-identification, 
the line thickness of the orbits is the same as for the portion 
of the characteristic curves in Fig. ^ to which they belong. 
Orbits in the lower and upper halfs of the diagram give rise 
to respectively stable and unstable orbits at low eccentricities 
when the bar is added. 



within bars, only few papers (Athanassoula et al. |1983 ; 
Contopoulos & Grosbol|l98|; Sellwood & Wilkinson p93 ) 
discuss them in discs surrounding bars. 

A good approach to investigate periodic orbits in a 2D 
rotating barred potential is to start with the axisymmetric 
limit. In this case, the only orbits which close whatever the 
value of rip are the circular orbits. All other bound orbits 
can be thought as a libration motion around these orbits 
and look like a rosette which never closes, except for some 
exceptional potentials like a point mass with Op = 0, or 
at resonances, where the radial and azimuthal frequencies 
ajR and uj^ satisfy the relation: 



(8) 



with integer values of nn > and n^. In the rotating 
frame, an nn/n^ resonant orbit closes after nn radial oscil- 
lations and 1 71,0 1 orbital periods. Outside corotation, is 
negative and one may speak of outer n/j/jn^l resonances. 
This paper discusses only outer resonances and the mi- 
nus sign in their labelling will be omitted. While in the 
axisymmetric case the orientation of the resonant orbits 
is arbitrary, the virtual introduction of a bar will retain 
only those orbits which are reflection symmetric with re- 
spect to (at least one of) the bar principal axes, say the 
X- and y-axes. Figure ^ displays several orbits of this kind 
with \n^\ = 1. There exists an infinity of resonant orbits 
which accumulate at corotation as nn — > oo. Those with 
ri/j > 6 or closing after more than one rotation will not 
be considered in this paper. For even n^, there exists two 
distinct resonant orbits for each value of the Hamiltonian, 
depending on whether the a;-axis coincides with orbital 
apocentre or pericentre. For odd n^j, there are twice as 
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Fig. 5. a) Characteristic diagram in the rotating axisymmetric potential (Eq. (^) for the circular orbit (circ.) and the lowest 
order outer nn/l resonant orbits which are symmetric with respect to the coordinate axes. The curves give the a;-coordinate 
normalised by the OLR radius of the orbits when they cross the a;-axis with y > as a function of the Hamiltonian value relative 
to i/oLR, the //-value of the circular orbit at the OLR, and normalised by v^^ assuming that the potential rotates clockwise in 
the inertial frame. The short-dashed line gives the zero velocity curve (ZVC) and the vertical dotted line the //-value of the 
orbits sketched in Fig. ^ The grey 3/1 and 5/1 resonance curves are shown only in the lower frame, which is a magnification of 
the rectangular box in the upper frame, b) Corresponding characteristic diagram in the barred potential of Eq. (fel) with a bar 
strength F = 0.10. The major axis of the bar coincides with the y-axis. The full and dotted parts of the curves stand for stable 
and unstable orbits respectively. The orbit labels refer to the nomenclature of Contopoulos & Grosbol (1989). The symmetric 
and asymmetric 3/1 and 5/1 resonance curves are plotted with grey lines in both the upper and lower frames, but labelled only 
in the lower one. c) Same as former diagram, but with F — 0.20. 



many solutions because the orbits are no longer reflection 
symmetric with respect to both axes. 

The characteristic curves of these resonant orbits in 
the H — X plane (Fig. ||a) all intersect the circular orbit 
curve (COC) at their point of lowest H, corresponding to 
a bifurcation. For even resonances, four branches emanate 
from the bifurcation, two from the COC and two from 
the resonance curve. The resonance branches above (to- 
wards larger x) and below the COC represent orbits with 
respectively apocentre and pericentre on the x-axis. The 
lower branch always passes through the zero velocity curve 
(ZVC), where the orbit becomes cuspy on the a;-axis and 
then develops loops at higher value of the Hamiltonian. 
For such loop orbits, the x-coordinate of the character- 
istic curves does not trace the pericentre but the place 
where the orbit self-intersect and x ^ Q. For odd reso- 
nances, the bifurcation has six branches: the two from the 
COC, two for the resonant orbits with radial extrema on 



the X-axis and which have properties similar to the former 
even resonance branches, and two for the resonant orbits 
with those extrema on the y-axis. The two latter branches 
have opposite x but degenerate into the same curve in the 
H — X characteristic diagram. All periodic orbits in the 
symmetry plane of an axisymmetric potential are stable. 

Figures and |^c show how the characteristic curves 
are modified when the bar component with major axis on 
the y-axis is added to the potential. The changes mainly 
occur at the bifurcations of the axisymmetric case. The 
bifurcations of the even resonances become gaps, with the 
right (low H) COC branch deviating into the lower reso- 
nance branch, and the upper resonance branch into the left 
COC branch, giving rise to a sequence of continuous orbit 
families. In the terminology introduced by Contopoulos & 
Grosbol ( 1989 ), the outermost of these families is called 
a;i(l), and the other families are divided into an upper 
x\(i) and a lower xi{i) sub-family at or near the point of 
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x,{2) 




drawn with thick Unes are examples of the perpendicu- 
larly oriented orbits that replace the circular orbit near 
the OLR radius and which have been proposed as a pos- 
sible explanation for either the Hyades and Sirius streams 



1 



Fig. 6. Orbits of the (left) and xi{2) (right) families in 

the model with bar strength F — 0.10, with the thick line repre- 
senting the orbit of lowest apocentre and of lowest Hamiltonian 
respectively. The shaded ellipse sketches the bar. The full and 
dotted lines represent stable and unstable orbits respectively. 
The nearly circular 2;i(l) orbits extend out to infinity. Note 
the unstable a;i(l) orbit at the 1/1 resonance. 



minimum H , where the stability of the orbits appears to 
reverse. The six-branch bifurcations of the odd resonances 
(see the 1/1 and 3/1 resonances in the figures) split into 
two pitchfork bifurcations, one involving the resonant or- 
bits symmetric with respect to the bar minor axis, which 
are stable near the bifurcation, and the other the resonant 
orbits non-symmetric relative to this axis, which are un- 
stable near the bifurcation and qualified as asymmetric. 
As a by-product of this splitting, the segment of the xi{i) 
characteristic curve between the two new bifurcations be- 
comes unstable. 

At low eccentricity (i.e. small H), only those orbits 
with a pericentre on the bar minor axis seem to be stable. 
This could be the influence of the Lagrangian points L4/5, 
which are stable fix points in the models. At high eccen- 
tricity, all orbit families undergo a change of stability and 
the characteristic curves in Fig. || are all interrupted at 
the point where the minor axis loop of the orbits touches 
the centre R = 0. The curves actually go beyond this 
point, but the resonance number changes. For instance, 
the 2/1 orbits become 2/3 orbits. When increasing the 
bar strength (from Fig. ||b to Fig. |^), the resonance gaps 
between successive xi{i) curves and the separation be- 
tween the pitchfork bifurcation pairs increases. The char- 
acteristic curves also become more twisted and the frac- 
tion of stable orbits decreases. In particular, almost all 
a;i(3) orbits are unstable for F = 0.20. It should be noted 



that other authors (Athanassoula et al. |1983| ; Sellwood & 
Wilkinson 1993| ) report more complicated characteristic 
curves for the Xi{2) and higher order orbit families near 
the ZVC in their models, probably as a consequence of an 
m — A Fourier component in the potential. 

Some orbits of the above families are plotted in D2000, 
but obviously missing all eccentric a;i(l) and xi(2) orbits 
with loops on the bar minor axis. A more complete set of 
orbits from these two families are given in Fig. |^. These 
orbits are the most important even ur periodic orbits be- 
cause they are usually associated with the largest invari- 
ant curve islands in surface of section maps. The orbits 



(Kalnajs |199lD or the Hercules-LSR bimodality (D2000) 
in the observed u — v distribution. As we shall see in the 
next sections, the stable eccentric xi(l) orbits also play an 
important role in shaping the local velocity distribution. 
Regular orbits trapped around them arc indeed unlikely 
to be heavily populated by stars, but represent forbidden 
phase space regions for chaotic orbits. 

The space coverage of the orbit families can be deter- 
mined from X — y plots like in Fig. ^. For each position 
in real space, there will be an (infinite) set of periodic or- 
bits passing through. The velocity trace in planar velocity 
space of the above described orbits, as well as the curves 
delineating some of the main resonances in the underlying 
axisymmetric potential^, are indicated in Figs. |^ and ^for 
various azimuthal angles ip and a realistic range of galac- 
tocentric distances relative to the OLR, and for two differ- 
ent bar strengths. Here the angle ip is measured from the 
bar major axis and increases towards the direction oppo- 
site to galactic rotation, i.e. coincides with the traditional 
in-plane inclination angle of the bar. All space positions 
are reached by many of the considered orbit families, and 
sometimes several traces are produced by orbits from the 
same family: for instance, at ip = 90°, there is a large range 
of R with three traces from a;i(l) orbits, mainly due to the 
loops on the a;-axis of the high-_ff orbits. The traces of the 
1/1, 1/1 asym, a;i(l), a;^(2), xi(2) and a;|(3) orbits with 
non-nearly vanishing |u|-velocity all fall very close to the 
associated resonance curves, indicating that the axisym- 
metric approximation used to compute these curves works 
well for nfl < 4. Not all resonant periodic orbits, i.e. those 
with traces on the resonance curves, are unstable. In par- 
ticular, orbits on the 2/1 (OLR) and 1/1 resonance curves 
are stable xi{l) and 1/1 orbits at 1^ = 90° and unstable 
a;*(2) and 1/1 asym orbits at = 0°, except for a xl{2) 
and a a;i(l) orbit with u = Q for R/Rqi^b. > l-l- Hence res- 
onance regions of phase space are not necessary unstable 
and depleted as asserted in D2000 for the OLR. 

There are sometimes several orbits from the same fam- 
ily plotted very close to each other, like for example the 
three unstable xi(2) orbits at -R/i?oLR = 0.9, ip = 30° and 
(v/vo,u/vo) ~ (-0.85, -0.80) for F ^ 0.20. This happens 
when the sequence of orbits within the family reverses its 
progression in the x — y plane towards a given direction 
and very close to the current space position, causing an ac- 
cumulation of orbits near this position with different local 
velocities. One may also note that the continuous transi- 
tions between some orbit families can cause periodic orbits 
from different families to have almost identical traces, as 
for example the x\{2) and xi(2) near the centre of the 
frame i?/i?oLH, = 1.1 and Lp = 90° in Fig. ||. 



2 As pointed out in D2000, the figures 2 and 4 in Fux ( |2000| ) 
wrongly display the OLR as a contour of constant H. This 
mistake is rectified here throughout the paper. 
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Fig. 7. Liapunov divergence timescale of the orbits in the u — v plane as a function of position in real space, for a bar strength 
F = 0.10. The horizontal and vertical axes of each frame are v/vo and u/vo respectively, with v positive in the direction of 
rotation and u towards the anti-centre, and the origin at the circular orbit of the axisymmetric part of the potential $o. The 
timcscalcs arc in units of local circular period tn in $0 and grcyscalc coded. The dark and white regions respectively represent 
regular and chaotic orbits. The oval and polygonal symbols indicate the positions of the periodic orbits, with a different symbol 
for each orbit family. Pull and empty symbols respectively stand for stable and unstable orbits. The full lines open towards 
the right (increasing v) are the H12 and H45 contours. The dash-dotted, solid, dashed and dotted lines open towards the left 
respectively give the locations of the outer 1/1, 2/1, 4/1 and 6/1 resonant orbits in $o- 
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Fig. 8. As Fig. Q but for a bar strength F = 0.20. 



The rather large number of stable simple periodic or- 
bits through each space position and the numerous stellar 
streams observed in the Solar neighbourhood may sug- 
gest that at least some of the m are r elated to periodic or- 
bits, as anticipated in Kalnajs ( 1991 ). For instance, beside 
the idea of the and low-eccentricity xi{2) induced 

streams near the OLR radius, the ~ 30° frames in Fig. ^ 
betray interesting stable eccentric xi{2) and 5/1 asym or- 
bits at R/RoLR = 1-2 and R/Rohn. — 1-1 respectively and 
with v/vo « —0.6 and u/vo ~ 0.05 — 0.15, which fall very 
close to the velocity of the young Arcturus stream (see 
Table llj). 



6. Liapunov exponents 

The next step after the periodic orbit search is to deter- 
mine the phase space extent of the regular orbits trapped 
around the stable closed orbits and of the chaotic orbits, 
which have no other integral than the Jacobi integral. The 
Poincare surface of section method is well suited to high- 
light the regular and chaotic regions in phase space at 
constant value of the Hamiltonian, but not at constant 
position in real space. A better tool for this purpose are 
the Liapunov exponents, which also allow to quantify the 
degree of stochasticity of the orbits. These exponents de- 
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scribe the mean exponential rate of divergence of two tra- 
jectories initially close to each other in phase space and 
are defined as: 



A(a;o, Aa;o) 



linr im^. 

Axo -^0 t \AXo\ 

t OG 



(9) 



where Xo and A the initial (t — 0) position and de- 

viation, and Ax(t) is the deviation at time t. Such limit is 
proven to exist and is finite for all bound orbits (Oseledec 
1968| )^. The value of A generally depends on the direc- 
tion of the deviation Axo and, if N is the dimension of 



ph ase sp ace, one can show (Oseledec |1968| Benettin et 
al. 1976) that there exists in fact N discrete exponents 
Ai > A2 > ... > Xn- However, the largest exponent Ai 
is the most important because it results from almost all 



deviations in Aajo-space (e.g. Udry & Pfenniger 1988D , 
and therefore is the one found in practice when the ini- 
tial deviation is chosen randomly. Moreover, this exponent 
exclusively determines the orbital stability: if Ai = 0, all 
exponents vanish and the orbit is regular, and if Ai > 0, 
the orbit is chaotic and the amount of chaos increases with 
Al. 

The numerical computation of the Liapunov exponents 
faces some problems related to the limits in Eq. (|^). First, 
the finite initial deviation Axo may rapidly grow as large 
as the size of the orbits themselves, especially for chaotic 
orbits, and thus Ax{t) must be occasionally scaled down 
by a large factor. Noting Aa;°(i) the deviation before the 
first rescalin g, this is done in a way similar to Contopoulos 
& Barbanis (1989), by normalising Ax{t) to |Aa;o|: 



Ax\U) 
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(10) 



every time when I Aa;* ^|> 1.3-10 "^/Jolr, so that after 
n iterations: 



\Axit)\ \Ax-it)\ 
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(11) 



Another problem is that one cannot integrate orbits over 
an infinite time and thus the time limit in Eq. (^ must be 
replaced by the value at some finite time, which has been 
set to about 3 Hubble times when the models are scaled 
to realistic physical units. 

The Liapunov exponent Ai has been calculated on a 
Cartesian grid of planar velocities for different positions of 
the observer, using the 2D potential of Eq. (||) and with an 
initial deviation of -1-1.3 • lG^^^i?oLR hi the i?-coordinate 
(Figs. 0, ^ and [ic| ). The results are presented as a diver- 
gence timescale t\ = 1/Ai in units of local circular pe- 
riod in the axisymmetric part of the potential, which pro- 
vides a more obvious and useful measure of stochasticity. 

^ In a realistic disc galaxy, the chaotic orbits not confined 
within corotation, and in particular the hot chaotic orbits, are 
not really bound but their escaping timescale is much larger 
than the Liapunov divergence timescale (see also Sect. |^. This 
is especially true in the case of our infinite mass logarithmic 
potential. 



Diagrams have been computed for every 10° in azimuth 
for R/Ro-LK = 0.9, 1.0, 1.05, 1.1, 1.2 and also at V3 = 25° 
over a larger radial range, but those diagrams between our 
30° sampling and at -R/-Rolr = 1.05 will not be shown 
here (and the same is true for the velocity distributions 
of the next section). We shall refer to such diagrams as 
"Liapunov" diagrams. 

The first thing to notice in these diagrams is that the 
stable and unstable periodic orbits fall within regular and 
chaotic regions respectively, as expected. There are some 
apparent exceptions like the 1/1 asym orbit at (/? = and 
large R which must owe to the limited velocity resolu- 
tion of the diagrams {— 0.04wo; see Fux 2001b for some 
higher resolution Liapunov diagrams). The fraction of 
chaotic orbits also obviously increases with bar strength. 
Furthermore, as a consequence of the four-fold symmet- 
ric barred potential, the diagrams at = and 90° are 
symmetric with respect to u = and, more generally, dia- 
grams at supplementary angles are anti-symmetric to each 
other in u, i.e. t\{R, 180° — ip, v, u) = t\(R, ip, v, —u). 

Ai Lp — 90° and for the radial range explored in Figs. |^ 
and H, the 2/1 and 1/1 resonance curves, as well as all 
other not highlighted resonances of the form 2/|n0| with 
integer \n^\^ are embedded in the middle of broad regu- 
lar orbit arcs separated by mainly chaotic regions which 
come closer to w = as the bar strength increases. For 
R/RoLB. ^ 1.1 (and ip = 90°), the regularity of the OLR 
arc gets destroyed near u = 0, leaving the place to an 
unstable a;*(2) orbit. Between the dominant regular orbit 
arcs, secondary regular arcs associated with resonances of 
the form A/\n^\ with odd integer \n^\ can also be identi- 
fied, especially for F = 0.10. This includes in particular 
the 4/1 resonance arc visible for R/RohB. 1-0. In fact, 
the chaotic regions between the broad resonance arcs are 
probably densely filled by tiny arcs of higher order regular 
resonant orbits, but the filling factor must be very low. 

At <p — 0, the situation is reversed: the 2/|n0| reso- 
nance curves lie in chaotic regions at large H which are 
spaced by regular orbit arcs right between the resonances. 
At intermediate angles, the regular and chaotic regions be- 
come offset from the resonance curves and the u-symmetry 
breaks. In particular, for (p ~ 30° and R/RohR ^ 1.0, i.e. 
realistic positions for the Sun, an extreme case of asym- 
metry arises near the OLR: a prominent region of regular 
orbits extends down to negative u, bounded roughly by the 
OLR curve on the right and penetrating well inside the hot 
orbit zone, whereas the positive u part of the OLR curve is 
surrounded by a wide chaotic region extending somewhat 
inside the H12 contour and coinciding very well with the 
u — V location of the Hercules stream. 

Figure |9| shows the w — u region occupied by the regular 
orbits trapped around the stable a;i(l) and xi(2) periodic 
orbits. To construct this figure, we have first derived many 
surfaces of section at mainly constant Hamiltonian inter- 
val AH = 0.054^0 to locate the islands of invariant curves 
around these periodic orbits. Then the orbits within each 
islands of these maps have been sampled by 50 regularly 
spaced points along a straight line across the central pe- 
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Fig. 9. Location in the u — v plane of the regular orbits trapped around the stable a::i(l) periodic orbits for a bar strength 
F = 0.20 (dark points) and around the stable Xi{2) orbits for _F = 0.10 (grey points). The computational details are described 
in the text. The horizontal and vertical axes of the frames and the different curves are as in Fig. ^. The l/-contours are given 
for F = 0.10. 



riodic orbit and within the outermost invariant curve of 
the island. Finally, all the resulting initial conditions are 
integrated for 20 rotations in the rotating frame and the 
velocities are plotted when the orbits pass within a dis- 
tance of dmax = i?o/80 of the considered space positions. 
The striped appearance of the regular a;i(l) and xi{2) re- 
gions in the diagrams are due to the discrete 7?-sampling 
and the broadening of the stripes to the finite value of 
do (not to an inaccurate orbital integration). The islands 
in the surfaces of section generally contain sub-resonances 
which have been included. The a^i(l) islands however have 
been truncated at the 1/1 resonance, so that the part of 
the xi{l) region on the high-w side of the 1/1 resonance 
curve is not represented in the u — v diagrams. The xi(l) 



regions are derived for F — 0.20 because the high level of 
chaos at this bar strength makes it easier to distinguish 
the boundaries of these regions in the surfaces of section, 
and the Xi{2) regions for F = 0.10 in order to emphasise 
the more regular case where these regions are larger. 

From this figure, it is obvious that the regular orbit 
arcs near the OLR, and in particular the prominent regu- 
lar region at (/? ~ 30° and R/Rohn ~ 1-0 discussed previ- 
ously, are produced by the regular orbits around the stable 
periodic orbits of the a;i(l) family, with an eccentricity 
increasing towards larger Hamiltonian values. The regu- 
lar regions associated with the stable xi(2) orbits can be 
viewed as divided into two parts, one involving only low- 
eccentricity orbits and the other one the higher eccentric- 
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ity orbits falling close to the 4/1 resonance curve. The low- 
eccentricity orbit part exists for R/Rqi^b. ^ 1-1 and over 
an angle range around ip = 90° increasing as R decreases, 
and is generally enclosed between the Hi2 contour and the 
OLR curve. Elsewhere, it is dissolved and only an unsta- 
ble x\{2) orbit remains (Fig. 0). The higher eccentricity 
part connects the low-eccentricity part at R/Rqlr ~ 0.9 
and then progressively detach from it as the 4/1 resonance 
curve moves away from the H12 contour at larger R. In 
particular, for i?/i?oLH, ^ 1-1 and for F = 0.10, there 
are no regular quasi-xi(2) orbits between the H12 contour 
and azimuthal velocities less than v/vo = —0.35, what- 
ever the angle (/?, and no low-eccentricity such orbits at all 
for (fi <i 50° . Hence for this bar strength, the Hercules-like 
mode found in D2000's simulations at R/Rqlr ~ l.lQand 
if 25 cannot be related to such quasi-a;i(2) orbits. 

Figure |l0| provides Liapunov diagrams over a larger 
radial range at 1^ = 25° and for F = 0.10. These dia- 
grams clearly show that the H12 contour marks the aver- 
age transition between regular and chaotic motion. Hence 
disc orbits are essentially regular and hot orbits essentially 
chaotic. The diagrams in Fig. |l^ also nicely illustrate the 
dependence of the H12 and H4C, contours with galacto- 
centric distance discussed in Sect. |[ and show how the 
velocity spacing between these contours decreases with in- 
creasing R. 



Martinet & Raboud (199!;) have computed a diagram 



More precisely at Rolh/R ~ 0.9. 



A(v,u) representing the relative pericentre deviation be- 
tween planar orbits integrated in a barred TV-body model 
and the corresponding orbits in the underlying axisym- 
metric potential, starting from a realistic space position of 
the Sun. Their diagram correlates well with our Liapunov 
diagrams in the sense that the larger values of A coincide 
with shorter Liapunov timescales. In particular, a clear 
jump of A occurs at H fa H12, with the high-i? side dis- 
playing much larger pericentre deviations on the average, 
and there is also a tail of small A-values at negative u 
extending inside the hot orbit region. 

It is worth mentioning that in a two-dimensional ax- 
isymmetric potential, all orbits are regular, i.e. have van- 
ishing Ai . Hence the chaotic regions discussed here are all 
produced by the influence of the bar alone. Also, the di- 
vergence timescale in the chaotic regions may be as low 
as a few orbital times. This property may have important 
consequences on the early evolution of the distribution 
function in barred galaxies, as we shall see in the follow- 
ing section. 



7. Phase space crowding 

Given the orbital structure of phase space, we now want 
to know how nature populates the available orbits. This 
is done resorting to test particle simulations with the fol- 
lowing two integral initial distribution function (Dehnen 
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1999a) 



fo{E,L,_) = 



■ exp 



(12) 



k{Rc) 2'Kd{Rc 

where E and Lz are respectively the energy and the z- 
component of the angular momentum, Rc{E) and Lc{E) 
the radius and angular momentum of the circular or- 
bit with energy E, and k the circular and epicycle 
frequencies, and and (y{R) the approximate sur- 

face density and radial velocity dispersion profiles. This 



is a modified Shu (1969) distribution where the radius 
of the guiding centre Rc is deduced from the energy 
instead of the angular momentum, with the main ad- 
vantage that the density function extends smoothly to- 
wards negative Lz- Adopting S(i?) ~ exp {~R/hu) and 
(T^(i?) = (Tq exp (— (i? — i?o)//i<T), where Ro is the galac- 
tocentric distance of the Sun and CTq the approximate ve- 
locity dispersion at that distance, the initial distribution 
function has three free parameters, which unless other- 
wise specified are set to hji/Ro = 0.33, (Tq/vo = 0.2 and 
h„/Ro = 1.0. This results in exactly the same initial con- 
ditions as for the default flat rotation curve case in D2000. 
In the beginning of the simulations, the non-axisymmctric 
part of the potential is gradually switched on from no 
contribution at t = to its full value aX t — 2th, where 
th = 27r/ilp is the rotation period of the bar, exactly the 
same way as in D2000 for the simulations with default 
integration time. 

For the time integration, D2000 uses a subtle back- 
ward integration technique based on the conservation of 
the phase space density along the orbits in coUisionless 
systems. The idea is to integrate back in time until t = 
the phase space points on a Cartesian grid of u — v ve- 
locities at a given space position {Ro,ip) and time tend, 
to determine the energy E and angular momentum Lz of 
the originating orbits in the initial axisymmetric poten- 
tial, and from these infer /(tend, ^o, </J, w) = fo{E,Lz). 
The advantages of this technique is that only the orbits 
strictly necessary to derive the evolved local velocity dis- 
tributions need to be computed, and that these velocity 
distributions are not affected by Poisson noise. Moreover, 
a unique simulation suffices to get u~ v distributions for 
different initial conditions, because /o comes in only after 
the orbit integration. 

Unfortunately, the backward integration technique 
faces two major problems illustrated in Fig. which 
shows the long term evolution of the planar velocity dis- 
tribution at Ro/RohR — 1-1 and — 25° using this tech- 
nique. The integration time in D2000 ranges from 4 bar 
rotation for most simulations, corresponding to only ~ 2 
orbital periods a,i R — Ro in the inertial frame, up to 
8 bar rotation for some cases, but always matching the 
growth time of the bar to half the total integration time. 
The frame at t = 4tb is similar to his results, reveal- 
ing a clear bimodal distribution. However, aX t = lOib, 
the valley between the two modes becomes heavily pop- 
ulated, destroying the bimodality, and at later times, in- 
curved waves appear in this valley with a spacing between 
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Fig. 11. Time evolution of the velocity distribution in the 
u — v plane at Ro/RohR — 1.1 and (f = 25° using the backward 
integration technique. The bar strength is F = 0.10 and the 
time is given in units of the bar rotation period th- 



the maxima decreasing with time. This is a typical signa- 
ture of ongoing phase mixing in a regular region of phase 
space, which here corresponds to the eccentric orbit part 
of the quasi- a;i(l) region according to the previous sec- 
tion. A very similar phenomenon, with similarly incurved 
waves, can also be barely recognised near the 1/1 reso- 
nance. Between the 2/1 and 1/1 resonances, phase mix- 
ing operates on a shorter timescale and the orientation 
of the wave fronts seems to change from nearly-vertical 
a.t t — 30tb to nearly-horizontal at t = 120tb- The back- 
ward integration technique in fact yields the fine-grained 
distribution function, which never smoothes out on suffi- 
ciently small scales, whereas the physical one to compare 
with the observations is the coarse-grained distribution. 
The second problem is related to chaos: at t ^ 30tb, the 
u — v distribution becomes noisy in the chaotic regions be- 
cause the phase space points integrated backwards from 
these regions sample the initial distribution function more 
randomly. Hence much longer integration times than in 
D2000 are required to obtain quasi-equilibrium (coarse- 
grained) distribution functions and to properly take into 
account the effect of chaos (remember that the divergence 
timescales for chaotic orbits is of the order of several or- 
bital periods), and one cannot escape the fate of smooth- 
ing the fine-grained distribution. 

Therefore, most of the test particle simulations in this 
paper were done by simple forward integration. The initial 
phase space density is sampled by N particles which are 
then integrated forward in time, and the u — v diagrams 
at space position {Ro , ^p) are constructed from all the par- 
ticles within a distance dmax = Ro/80 from that position 
(corresponding to dmax — 100 pc for Ro — 8 kpc), using 
a Cartesian velocity binning with a bin size of O.OOSuo- 
To increase the particle statistics, the u — v distributions 
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Fig. 12. Velocity distribution in the u — v plane as a function of space position in the test particle simulations with a bar 
strength F — 0.10. The distributions are time averages within the interval from 25 to 35 bar rotations after the beginning of 
the simulations. The velocity contours are as in Fig. hi whereas the axis labels, the Hamiltonian contours, the resonance curves 
and the periodic orbits are as in Fig. ^. 



are averaged over 10 bar rotations and then smoothed 
within a square of 11 x 11 bins. The time average, vifhich 
is hardly possible in the backward integration technique, is 
also very convenient to reduce the phase mixing problem, 
as the contribution of each particle to a u — v distribu- 
tion becomes proportional to the time the particle spends 
within the volume where the distribution is computed. 
With the forward integration technique, the time evolu- 
tion of the velocity distribution can be followed within a 
unique simulation. The test particle simulations of this pa- 
per all have A'' = 10^ and the velocity distributions have 
been derived within two time intervals, 25 < t/t^, < 35 
and 55 < t/t^j < 65. Since the distance parameters in /o 



scale as Ro and not i?oLR, the results at different Ro/ Rqlr 
require distinct simulations. 

Figure 12 shows the u—v distributions at various space 



positions averaged over the time interval 25 < t/t^, < 35 
for a bar strength F = 0.10. As for the Liapunov diagrams, 
the distributions are obviously symmetric with respect to 
u = for (p ~ and (p = 90° and the distributions at 
same radius but supplementary angles are anti-symmetric 
to each other in u. This is clearly not the case in the sim- 
ulations of D2000 (see his figure 2, where F — 0.089), pro- 
viding a further argument that these have not achieved a 
quasi-stationary regime. Moreover, the traces of the stable 
xi (1) periodic orbits away from the 2/1 resonance curve lie 
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Fig. 13. Same as Fig. |l^, but for a bar strength F = 0.20 and a time average over the interval 55 < t/th < 65. 



closer to the high angular momentum peak of the velocity 
distributions. For this bar strength and the adopted values 
of the parameters in the initial distribution function, there 
is also no clear bimodality with a deep separation valley 
at i?o/i?oLR ^ 1-1 and near = 30°, although a clear 
density excess remains at low v and positive i^. However, 
all space positions where a low-eccentricity regular Xi(2) 
region exists (see former section and Fig. ^ present a nice 
bimodality, with the low angular momentum mode coin- 
ciding very well with that region and always peaking inside 
the Hi2 contour. 

While the traces of the non-resonant xi{l) and xi{2) 
orbits are generally embedded within their associated 
quasi-periodic orbit modes, they do not necessarily ex- 

^ Doubling the average integration time can reinforce the 
bimodality, as shown in Fig. bda. 



actly coincide with the peak of these modes, especially in 
the xi(2) case (e.g. Ro/Rolr = 1-0 and ,^ = 60° - 120° in 
Fig. |l^). Furthermore, since the quasi-periodic orbits cover 
a larger space area than the periodic orbits themselves, 
quasi- a;i(l) and quasi-a;i(2) modes may occur even at po- 
sitions where no a;i(l) or a;i(2) orbit are passing through. 

Increasing the bar strength (Fig. ^3|) provides a bet- 
ter understanding of how the velocity distributions are 
affected by chaos. Now, there appears to be an obvious 
second source producing a low angular momentum mode, 
which adds to the quasi-a;i(2) orbit flow at the space re- 
gions reached by these orbits, and acts alone elsewhere, 
as for instance at Ro/Rolr = 1.1 and (p = 30°. The over- 
density in velocity space generated by this second source 
correlates very well with highly stochastic regions in the 
Liapunov diagrams (see Fig. ^) and seems to always peak 



16 



R. Fux: Order and chaos in the local disc stellar kinematics induced by the Galactic bar 



outside the H12 contour. At (p — 90°, the overdensity cul- 
minates at It « and is enclosed by the regular arc of 
eccentric quasi- orbits. At ip ~ 0, these quasi-a;i(l) 
orbits occupy the u — region and chaotic overdensities 
happen symmetrically on both positive and negative u 
sides of this region. At ip 30°, the quasi- a;i(l) region is 
located at negative u and there is a large chaotic overden- 
sity at positive u. 

These properties result from the decoupling between 
the regular and the chaotic regions of phase space. Since 
chaotic orbits cannot visit the regions of regular mo- 
tion and, vice versa, regular orbits avoid the chaotic re- 
gions, the distribution function in each of these regions 
evolves in a completely independent way. In the regular 
regions, it recovers roughly the initial distribution after 
phase mixing, whereas in the chaotic regions, it is sub- 
stantially modified through a process known as chaotic 
mixing and which operates on the Liapunov divergence 
timescale (e.g. Kandrup |200l| ): the particles on chaotic 
orbits quickly disperse within the easily accessible phase 
space regions, i.e. not impeded by cantori or an Arnold 
web, and converge towards a uniform population of these 
regions. The dominant manifestation of chaotic mixing 
is a net migration of particles from the inner to the 
outer space regions. For instance, in the simulations with 
Ro/RohR = l-li the scale length of the radial particle dis- 
tribution, which remains very close to an exponential in 
the range 0.5 ^ R/Ro ^ 1.25, increases by ^ 30% for 
F = 0.10 and by ~ 90% for F = 0.20 within this radial 
range. This migration is particularly marked for particles 
on hot chaotic orbits because the region inside corota- 
tion initially represents a large reservoir of such particles 
and because these particles can freely pass over corota- 
tion. As a consequence, in the explored Ro/RoLn range, 
the chaotic regions in the u — v diagrams are more heavily 
crowded than the regular regions a.t H ^ H12, therefore 
corresponding to local overdensities. 

At first glance, there seems to be a rather continuous 
transition from the quasi-xi(2) orbit mode to the main 
chaotic orbit mode when moving across the OLR radius 
towards increasing Ro , with always a single effective peak 
showing up and with the involved quasi-a;i(2) region pro- 
gressively dissolving in the chaotic one. But in some cases, 
the two mode-generating sources really contribute to dis- 
tinct peaks in the velocity distribution (see Fig. |20| d for 
an example). 

The process of chaotic mixing leads to velocity dis- 
tribution contours which are parallel to the contours of 
constant H in the chaotic regions (e.g. Fig. |l^). This 
property is also consistent with Jean's theorem stating 
that the distribution function in a steady-state system de- 
pends only on the integral of motions. The only integral 
for the chaotic orbits in the present 2D barred models is 
the value of the Hamiltonian, hence the distribution func- 
tion and therefore the corresponding velocity distributions 
at fixed space position should be a function of only H in 
the chaotic regions. It should be noted that the Jeans the- 
orem does not strictly apply to the hot and disc chaotic 



orbits. Indeed, these orbits are not energetically bound (in 
terms of H) and thus the phase space density around such 
orbits and within the finite space volume of the galaxy 
should decrease with time, conflicting with the steady- 
state assumption of the theorem. However, the escaping 
timescale of these chaotic orbits, which is essentially con- 
trolled by Arnold diffusion across the confining cantori, is 
much longer than the Liapunov divergence timescale and 
even the Hubble time, and thus the density in the phase 
space regions covered by these orbits can be considered as 



almost constant (see also Kaufmann & Contopoulos 1996 
and references therein). 

Secondary chaotic orbit overdensities also occur be- 
tween the a;i(l) and the 1/1 regular arcs, especially at 
Ro/RoLK > 1.0 and 30° < < 150° (Fig. These sec- 
ondary overdensities and the above described main over- 
densities connect to each other in phase space, i.e. are 
traced by the same orbits. Hence it is also expected that 
the u — v density at constant H-vahie is the same for all 
overdensities. This is only roughly the case in Fig. ^3[ 
probably because the smoothing of the diagrams lowers 
the peaks of the tiny secondary overdensities relative to 
the broader main overdensities. Small chaotic overdensi- 
ties may sometimes even be noticed beyond the 1/1 res- 
onance curve. However, at high angular momentum, the 
hot orbits spend most of their time in the outer galaxy, 
far away from the influence of the bar, and thus become 
more regular (the energy and angular momentum are more 
nearly conserved). The eccentric xi{2) regular regions can 
also represent density depressions between chaotic regions 
in the velocity distributions, as can be marginally inferred 
for exaniple from the Ro/Rohn = 0.9 and (p = 90° frame 
in Figs. I and [l2[ 

The valley between the main high-Lz velocity mode 
(or LSR mode after Dehnen) and the main chaotic orbit 
mode is generally close to the H12 contour, reflecting the 
decline of the hot orbit population a,s H ^ H12. Such 
a valley should in principle also exist between the LSR 
mode and the secondary chaotic overdensities (see for in- 
stance Ro/Rqlk = 1.1 and (p = 30° in Fig. |l^. The main 
chaotic orbit mode also seems to always peak between the 
H12 and H45 contours in Fig. |l^, but this is not true for 
all our test particle simulations, as demonstrated by the 
hn/Ro = 0.2 frame in Fig. and Fig. pO|a. However, this 
property might be more generic for self-consistent models 
(see Sect. 0). For some not fully understood reasons, the 
symmetry properties mentioned previously for the case 
F = 0.10 are somewhat less evident for F ~ 0.20, despite 
the longer integration time. 

D2000 attributes the valley between the main LSR 
mode and the Hercules-like mode to stars scattered off 
the OLR, in the sense that the resonance generates chaotic 
orbits. In particular, he claims that the unstable xl (2) or- 
bit falls exactly between the two modes. This is not quite 
correct for a stochastically induced Hercules-like mode, 
as is best illustrated by Fig. ^ for Ro/RohR ^ 1-0 and 
p ^ 30°, the M > part of the 2/1 resonance curve passes 
through the Hercules-like mode and the x* (2) orbit clearly 
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lies within the mode at Ro/RohR = 1-2, and the low den- 
sity region below this mode is due to regular resonant or- 
bits. D2000 also claims that in his simulations the exten- 
sion of the LSR mode to m < at u « — O.lvo is caused 
by the elongation of the (presumably quasi- a:i(l)) orbits 
near the OLR. Our results indicate that at least the fi- 
nal part of this extension, corresponding to the secondary 
overdensities, is produced by chaotic orbits. Such an ex- 
tension exits in the observations, but is only significant 
down to heliocentric u sa —60 kms^^. 

The velocity distributions at the two different mean 
integration times <t>~ 30ib and 60tb reveal some sec- 
ular evolution. As < t > increases, the crowding contrast 
between the regular and chaotic regions becomes more evi- 
dent, with denser high-iJ chaotic regions and deeper regu- 
lar region valleys. The quasi- a;i(l) mode squeezes towards 
its high angular momentum side for Ro/Rqlr 1.1 and, 
especially in the case F = 0.10, the peak of the quasi-xi(2) 
mode moves closer to the H12 contour for Ro = Rolr, 
betraying a longer phase mixing timescale near this reso- 
nance^ 

In a real galaxy, the presence of mass concentrations 
like giant molecular clouds and of transient spiral arms 
will prevent the strict conservation of the Jacobi integral 
and cause the stars to diffuse from the regular regions 
to the chaotic regions of phase space and vice versa (see 
also Sect. |o|). The chaotic regions should be very efficient 
in heating galactic discs and the communication between 
the two kind of regions may even allow to heat regular 
regions. The quantification of this phenomenon might be 
an interesting problem to study, but is beyond the scope 
of this paper. 

Finally, Figs. |lj and |l^ also suggest that with increas- 
ing bar strength, the velocity dispersion ratio ay/(Tu de- 
creases and, as reported by D2000, the u-range and in 
particular the mean u-velocity of the main chaotic over- 
density become larger. 

8. Changing the initial conditions 

How can the initial conditions be changed in order to in- 
crease the population of the hot chaotic orbits, and in 
particular the density in the Hercules-like chaotic over- 
density ? Figure ^ shows how the three parameters of 
the initial axisymmetric distribution function (Eq. (p^)) 
individually affect the final velocity distribution for a real- 
istic position of the Sun and F = 0.10, using exceptionally 
the backward integration technique which, as mentioned 
in Sect. |^, is very convenient for this purpose. A long time 
integration is adopted to reduce the phase mixing problem 

® This could be a consequence of the fact that the linear 
radial oscillation frequency of the quasi-periodic orbits around 
the least eccentric stable xi{2) orbits near the OLR radius, i.e. 
the non-axisymmetric analogue of the epicycle frequency for 
these periodic orbits, is close to the radial frequency wr of the 
Xi{2) orbits themselves. A similar argument may also hold for 
the slow phase mixing noticed in Fig. |l^ within the resonant 
part of the quasi- a;i(l) region. 
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Fig. 14. Velocity distribution in the u — v plane after 120 
bar rotations as a function of initial conditions and using 
the backward integration technique. The space position is 
-R0/-R0LR = 1.1 and (p = 25° and the bar strength F = 0.10. 
The left frame shows the result for the default values of the 
parameters, i.e. ao/vo = 0.2, hc/Ro = 1 and hn/Ro = 0.33, 
and the other frames the results when changing only one pa- 
rameter at a time to the indicated value. The velocity contours 
are as in Fig. ^ and the circular arcs represent the H12 and -ff45 
contours. 

and the u — v distributions are smoothed the same way as 
the previous ones based on the direct integration method. 
A comparison of the default parameter velocity distribu- 
tion with the corresponding distribution (at if = 30°) 
in Fig. |l^ indicates that both integration techniques give 
very similar results. 

Increasing the overall initial velocity dispersion (top 
frames in Fig. yields a larger final velocity disper- 
sion. Hence in this kind of simulations the particles re- 
member the initial conditions and the action of the bar 
is unable to completely erase them. Also, the u — v den- 
sity in the Hercules-like stream is enhanced relative to the 
density within the main velocity mode. This is because 
the larger velocity dispersion lowers the peak of the lat- 
ter mode, and because a larger CTo increases the average 
Hamiltonian value of the particles and thus the fraction of 
hot chaotic particles. Reducing the initial velocity disper- 
sion scale length keeping the same velocity dispersion at 
Ro (middle frames) renders the inner regions hotter and 
hence mainly increases the density of the velocity distri- 
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Fig. 15. For two different space positions and after 120 bar 
rotations, traces in the u — v plane of the trajectories which 
were in OLR with the unperturbed rotating potential at f = 
(bar strength F = 0.10). The points are colour-coded according 
to the angle a between the line joining the apocentres of the 
initial resonant orbit and the major axis of the vanishing bar 
potential: black if |q:| < 45°, light grey if |q| > 70°, and dark 
grey for intermediate angles. 



button at low angular momentum. To increase the relative 
fraction of stars in the Hercules-like stream, the most ef- 
ficient way seems to start with a smaller disc scale length 
(bottom frames). This represents a higher initial space 
density of the disc in the inner regions where the parti- 
cles have larger values on the average and thus again a 
larger fraction of hot chaotic particles which will be spread 
over the whole disc by the barred potential. 

By changing the initial conditions, it is therefore pos- 
sible to achieve a more pronounced main chaotic overden- 
sity than in the results based on the default parameters 
and also to match more precisely the observed velocity 
dispersions in the Solar neighbourhood. 

9. Resonances and stochasticity 

A worthful exercise is to determine what happens to the 
periodic orbits which are initially in (outer) 2/1 reso- 
nance before the growth of the bar. This is illustrated in 
Fig. ^ which highlights for two different space positions 
at t = 120tb the points in the u~v plane corresponding to 
trajectories which were on such orbits at t = 0. The dia- 
grams are built by integrating backwards the trajectories 
passing through a Cartesian u — v grid and marking all 
the points on this grid originating from initial orbits with 
{ui^ + 2uifi — r2p)/r2p < lO^'^. The darkness of the points 
reflects the angle a between the major axis of the initial 
resonant orbit and the major axis of the then vanishing 
bar potential, with darker points standing for smaller an- 
gles, i.e. resonant orbits with apocentre closer to the bar 
major axis. 

At = 90° and i?/i?oLR — 0.9, there is a wide region 
of regular orbits around most of the 2/1 resonance curve 
(see Fig. 0). The trajectories with small a's clearly fall 



in the inner part of this region, while those with larger 
a's appear rather on its edge and are spread within the 
neighbouring chaotic regions. In addition to the regular 
versus chaotic phase space decoupling, the fact that the 
OLR curve is associated with a valley in the planar ve- 
locity distribution (see Fig. |l2|) is also because the phase 
space density in the chaotic regions bounding the regular 
orbit arc around this curve is forced to be a function of 
H only, hence lowering the density on the low-u side and 
increasing it on the other side relative to the initial den- 
sities. At = 25° and i?/i?oLR = 1-1, the part of the 2/1 
resonance curve within the u > chaotic region has no 
nearby small-a points and is embedded in a broad cloud 
of high-a points. 

Hence, the initial 2/1 resonant orbits more nearly 
aligned with the bar major axis end into regular regions, 
while only those more nearly perpendicular to this axis 
are scattered into chaotic regions. This is consistent with 
the stability properties of the low-eccentricity a;i(l) and 
a;* (2) orbits in the full barred potential, which are both 
2/1 periodic orbits. 

10. A^-body models 

In addition to the test particle simulations, we have also 
run a high-resolution A^-body simulation whose predic- 
tions can be compared with the results of the former sec- 
tions. The main difference of this simulation with respect 
to the previous simulations is that it is three-dimensional 
and completely self-consistent, i.e. with no rigid compo- 
nent and allowing the development of spiral arms. The de- 
scription of the simulation hereafter is based on physical 
units in which the initial conditions provide a reasonable 
axisymmetric model of the present Milky Way. In these 
units, i?oLR ~ 9 kpc dXt — 2.32 Gyr. 



The simulation, also discussed in Fux (200C ), is similar 
to the simulations presented in Fux ( 1997 ), starting from 
a bar unstable axisymmetric model including a nucleus- 
spheroid, a disc and a dark halo component with parame- 
ters set to a = 1 kpc, Mns = 2.6 • lO^" Mq, /i^^ = 3.2 kpc, 

= 300 pc, Md = 5.0 • 10^° Mq, h = 9.1 k pc an d 
A/dh = 2.6-10" Mq (in the same notation as in Fux |l997|) . 
It involves 14 299 552 particles of which 5 553 784 belong to 
the disc component. The simulation uses the double polar- 
cylindrical grid technique described in Fux ( 1999D to solve 
the gravitational forces, with Nr x A^^ x A^x — 94 x 96 x 253, 

— 25 pc and imposed reflection symmetry witli respect 
to the plane z = and the z-axis. The time integrator is a 
leap-frog with a time step At = 0.05 Myr. The simulation 
has been carried out until t — 2.65 Gyr. The phase space 
coordinates have been saved every 50 Myr for all particles 
and every Myr for the disc particles within a fixed annulus 
7.5 < i?/kpc < 10.5. 

After the formation of the bar at about 1.4 Gyr, the 
simulation reveals a complex velocity structure with mul- 
tiple streams occuring almost everywhere in the disc. 
The streams, as well as the velocity dispersions, remain 
very time dependent, even within space regions corotat- 
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Fig. 16. Some properties of the A'^-body simulation averaged 
over 2 < t/Gyr < 2.65 and scaled to a time independent -Rolr 
(see text for details), a) Face-on surface density of the lumi- 
nous (disc -|- nucleus-spheroid) particles, with contours spaced 
by 0.5 magnitude, b) Effective potential in the plane z = 0, 
with the spacing of the contours increasing by a factor 1.2 as 
one moves from L4/5 to lower "I>off regions. The curves made 
of large and small dots represent respectively the azimuthal 
minima and maxima of <l>cff at local radius. In these first two 
frames, galactic rotation is clockwise and the crosses indicate 
the Lagrangian points L1/2 and L4/5. c) Effective potential 
in the vertical plane passing through the Lagrangian points 
Li and L2 and shown by the dotted straight line in the for- 
mer frame, d) Rotation curve of the axisymmetrised average 
model, with uolr = i'c(-Rolr). e) Radial dependence of the 
azimuthal and radial bar strengths Ftp and Fr. The unlabelled 
dotted curve is the azimuthal strength of the analytical barred 
potential given by Eq. (Eh for F = 0.20. 



ing with the bar. Some mpeg movies showing the time 
evolution of the u — v distribution for a realistic posi- 
tion of the Sun are available on the web at the address 

However, 



http : //www.mso . anu. edu. au/'f ux/streams .html 



the strong time dependency is certainly partly the conse- 
quence of incomplete phase mixing as in the test parti- 
cle simulations of the previous sections. Therefore, we will 
again average in time the w distributions resulting from 
the A'^-body simulation to minimise this problem and thus 
focus on the average properties of these distributions. 

A complication arises from the fact that the pattern 
speed of the bar decreases with time due to the angular 
momentum the bar loses to the (live) dark halo, so that 
the OLR does not lie at a fixed radius anymore but moves 
outwards. The decrease of ilp is probably not that large 
in real disc galaxies with a substantial gas fraction like the 
Milky Way, as shown by self-consistent numerical s imula - 
tions with an SPH component (e.g. Friedli & Benz 1993| ). 
Hence, the u — v distributions at a given Ra/ Rolr must 
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Fig. 17. a) Characteristic diagram for the low-eccentricity 
xi{l), xl{2) and xi{2) periodic orbits in the time averaged 
A'^-body model. The full and dotted lines indicate stable and 
unstable orbits respectively, and the dashed line the zero ve- 
locity curve. voLR is the circular velocity at the OLR in the 
axisymmetrised potential. All curves are interrupted at high H 
(see text), and the a;i(l) curve is truncated on the low-H side 
at the outer 1/1 resonance, b) Sequence of a;i(l) orbits in the 
x~y plane for the same model. The orbits are sampled at con- 
stant Hamiltonian interval and cover the whole range of the 
corresponding characteristic curve in the former diagram. The 
orientation of the bar is sketched by the shaded ellipse and its 
rotation is clockwise. 

be computed within space volumes comoving with i?oLR- 
We found that the evolution of the bar position angle d{t) 
over the time interval 1.5 < t/Gyr < 2.65 can be well 
fitted by an analytical function of the form: 



exp(-^i), 



(13) 



where i^oo, ^^o and fi are the free parameters, with residu- 
als in I? never exceeding 5°. From this we obtain flp{t) = 
d-ff/dt and RcB_{t) = v'^/V,p{t), where v'^ is matched so 
that Rqyl — [^(^1/2) + -^(-^4/5)1/2 on the average. The 
radius of the OLR is then derived from i?cR via the flat ro- 
tation curve relation, which yields a very good approxima- 
tion of i?oLR on the average, and represents the smoothly 
changing reference adopted for the distance normalisation. 
Because the rotation curve remains nearly flat and con- 
stant with time at intermediate i?, no velocity scaling is 
required. 

The time average is done over the interval 2 < t/Gyr < 
2.65 and as underlying mass distribution and potential as- 
sociated with the average velocity distributions, we take 
those resulting from the sum of the 50 Myr spaced output 
models of the simulation within the same time interval and 
with the mass in each model rescaled as the distances to 
preserve the velocity scale. The number of models added 
together is rather low (only 14), yielding only a crude es- 
timate of the true average quantities, especially regarding 
azimuthal variations in the spiral arm regions. Figure ^ 
shows some properties of the resulting average model. At 
given radius _R, the azimuthal and radial bar strengths 
F^{R) and Fji{R) (Fig. [l^e) are defined as the most ex- 
treme values over azimuth of \A^{R,ip)/A'^^^^"^{R)\ and 
[An{R,^) - Af''y''\R)]/Af'''''''\R) respectively, where 
A(ji and An are the azimuthal and radial accelerations in 
the plane z = and ^^'"y™ jg ^j^g axisymmetric part of 
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Fig. 18. Time averaged planar velocity distribution of the disc particles in the A*'-body simulation as a function of position in 
space and within vertical cylinders of radius _Ro/80. The horizontal and vertical axes of the frames are v/vo and u/vo respectively, 
with Vo being the circular velocity at local radius in the corresponding axisymmetriesed average model. The velocity contours 
are as in Fig. as well as the H-contours (drawn at 2 = 0), the resonance curves and the traces of the periodic orbits, which 
all refer to the average rotating potential. Only orbits from the x\(l), 2^1(2) and x\[2) families and within the Hamiltonian 
range scanned by the characteristic curves in Fig. [l^ a are indicated. 



Ar- In particular, F^{a = i?cR/l-25) coincides with the 
definition of bar strength used in the former sections. The 
time averaged model has f(^( a) k. 0.2, i.e. a rather strong 
bar. However, Buta & Block ( ^001 ) have introduced a bar 
strength classification in terms of a parameter Qb, cor- 
responding here to the radial maximum of F^{R), and 
present galaxies with values of this parameter up to 0.65. 
Our average model has Qh ~ 0.325 and thus falls only in 
the middle of the range covered by real galaxies. 

The disc scale length between corotation and slightly 
outside the OLR (Fig. [l6^) has substantially increased 
with respect to the initial conditions, in agreement with 
the test particle results. The effective potential (Fig. 
indicates Lagrangian points which significantly lag the bar 
principal axes. This is the effect of the spiral arms starting 
at the end of the bar, which were absent in the test parti- 
cle simulations and now produce a twist of the po tentia l 
well. Note however that, as pointed out by Zhang ( 1996|) , 
there is a phase shift between the potential and the den- 
sity wells, with the former leading the latter outside the 
bar. The characteristic curves of the main periodic orbit 



families and the x — y configuration of the a;i(l) orbits in 
the average potential are given in Fig. The character- 
istic curves are truncated at large H , where they start to 
bend and interfere with each other presumably as a con- 
sequence of the sharp phase change in the potential well 
at large radii (see Fig. p^b). The a:i(l) orbits respond to 
the twisted potential well by having their major axis de- 
parting more and more from the y-axis as H decreases, 
and the closed orbits of the other families share a similar 
response. Since the cusps of the cusped orbits now occur 
away from the coordinate axes, the characteristic curves 
no longer reach the ZVC. 

The time averaged u — v distributions are presented 
in Fig. The diagrams are computed by summing the 
1 Myr spaced outputs of the simulation, yielding much 
more accurate time averages than for the properties based 
on the 50 Myr spaced outputs, and using the same space 
volumes, velocity binning and smoothing procedure as for 
the test particle simulations. However, contrary to the lat- 
ter simulations, the diagrams are based on a unique sim- 
ulation and therefore the initial conditions for each dia- 
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gram now scale as i?oLR instead of Ro . Hence the pattern 
speed and size of the bar are not the only parameters that 
change for different values of i?o/i?oLR- In particular, the 
initial velocity dispersions decrease with i?,, causing a simi- 
lar trend in the final velocity distributions. Diagrams have 
been derived at an azimuthal interval A(p — 10°, but only 
those at 30° spacing are shown. 

A priori, some properties inferred from the test parti- 
cle simulations seem less robust in the iV-body simulation: 
the non-resonant xi(l) orbits are somewhat less coincident 
with peaks in the velocity distributions, and the resonant 
a;i(l) orbits, i.e. those with traces on the 2/1 resonance 
curve, are less correlated with depleted u — v regions at 
high eccentricities. Moreover, the low angular momentum 
peak often lies well inside the hot orbit region even when it 
is still mostly associated with regular quasi-xi(2) orbits, 
as indicated by the presence of a stable low-eccentricity 
a;i(2) orbit. In fact, a closer inspection reveals that the 
velocity distributions in the iV-body simulation at given 
Ro/RoLB. appear to best match those of the test particle 
simulations at a typically 10% larger value of Ro/Rqlr- 
This can be explained by a delayed response of the phase 
space density distribution to the growing absolute OLR ra- 
dius: the u — v distributions do not instantaneously adjust 
to the moving i?oLR and therefore always reflect an earlier 
smaller i?oLR value instead of the current one. Hence the 
velocity distributions in Fig. |lj should virtually be shifted 
upwards by roughly one frame to be more consistent with 
the Ro/RoLR scale and the other plotted informations. 
Actually, the effective i?oLR must be a function of the or- 
bits. 

With such a correction in mind, the A^-body velocity 
distributions now share much more the same properties as 
in the test particle simulations. The low angular momen- 
tum mode, when purely induced by chaotic orbits as in 
the frame Ro/Rolr = 1.1 and (/j = 30° of Fig. |l|, also no 
longer peaks outside the H^^^ contour, but rather between 
the Hi2 and the i?45 contours. The velocity distributions 
most resemble those of the test particle simulations with 
a bar strength F = 0.15 (not shown in this paper, but 
see Fux 2001a ). Since F « 0.20 in our average A^-body 
model, this suggests that the velocity distributions are 
less responsive to the bar strength in the more realistic 
3D A^-body simulation than in the 2D test particle simu- 
lations. 

The symmetries found in Sect. |^ for the velocity dis- 
tributions, and in particular the it-symmetry at (f — 
and (p = 90°, obviously break and the velocity distribu- 
tions in the A^-body model at given ip and fixed effec- 
tive radius seem to compare best with the corresponding 
distributions in the test particle simulations at an angle 
~ (p—10°, suggesting that the velocity distributions know 
about the spiral arm induced average local twist of the 
potential well relative to the bar major axis. While this 
is especially true for the more regular \ow-H regions, the 
velocity distributions show no significant phase shift at all 
in the hot orbit region. The reason is because the hot or- 
bits are more eccentric and therefore are sensitive to more 
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Fig. 19. Comparison of the normalised Hamiltonian values of 
a random selection of 2 • 10^ disc component particles in the 
A-body simulation at times t = 2000 and 2500 Myr. volu and 
-ffoLR respectively are the velocity and Hamiltonian value of 
the circular orbit at the OLR of the axisymmetrised potential. 
The points in grey represent the particles inside the circle pass- 
ing through the Lagrangian points I/1/2 att — 2000 Myr. These 
particles exist down to and beyond the smallest displayed value 
of H, as the effective potential at the central Lagrangian point 
is (^effiLs) — -H'oLR)/t'OLR ~ —5.2. The solid lines indicate 
the normalised H12 and H45 values and the dash-dotted lines 
the normalised //-value of the circular orbit at the outer 1/1 
resonance (in the axisymmetrised potential). 



inner features of the potential. It should be noted that 
in A^-body simulations like the one presented here, spiral 
arms are particularly strong during about 1 Gyr after the 
formation of the bar, so that the reported effects of spiral 
arms are probably overestimates for the Milky Way if the 
Galactic bar is old. 

A potentially important difference of 3D models with 
respect to 2D models is that the effective potential a star 
experiences near corotation depends on its distance from 
the Galactic plane (see Fig. |l6|c). This raises the aver- 
age value of the Jacobi integral required for the stars to 
cross the corotation radius and thus renders such a cross- 
ing more difficult. For stars on the Solar circle, the higher 
effective value of H12 is not compensated by their depar- 
ture from the Galactic plane or a w velocity component. 
Indeed, in our average A^-body model, the change of effec- 
tive potential near corotation when moving from z = to 
z = 300 pc is A<i>off/wQLR ~ 0.233 (with Wqlr as defined 
in the caption of Fig. |l^ , while this change at the OLR 
of the axisymmetrised potential is only 0.004, and adding 
a vertical velocity of w/^olr = 0.1 only increases H/vq^^ 
by 0.005. Hence 2D models certainly exaggerate the traffic 
of stars on hot orbits travelling from one side of corotation 
to the other. 

Finally, Fig. |l9| shows how the value of the Hamiltonian 
is conserved during the A^-body simulation. The main re- 
sult is that the i/- values at different times are much better 
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related for bar particles than for (the dynamically defined) 
disc and hot particles. In particular, bar particles remain 
bar particles, but disc particles can easily transform into 
hot particles and vice versa for (-ff — ^^olr)/wolr ^ —0.2, 
supporting the presumption in Sect. ^ that spiral arms 
may induce exchanges between regular and chaotic phase 
space regions in real galaxies. Near the 1/1 resonance, the 
disc particles also reveal a smaller scatter of their past ver- 
sus present H relation. The normalised Hamiltonian (as 
well as the absolute non-rescaled Hamiltonian) increases 
on the average for the disc particles, which may be un- 
derstood by the fact that the contribution of the term 
— f2pi?^/2 to $cff diminishes as the bar rotates slower, and 
decreases for the bar particles, owing to the deepening of 
the central potential well. 

It would be interesting to investigate the evolution of 
the particle iJ-values in an A^-body simulation with a bar 
rotating at a constant frequency, for example without in- 
cluding a live dark halo component, in order to disentangle 
the effects of the spiral arms from the effects of a decreas- 
ing pattern speed. It would be also worth to explore the 
diffusion of particles from regular to chaotic phase space 
regions and vice versa, starting with 2D A^-body simula- 
tions in a first approach. One problem to be clarified is 
why the A^-body velocity distributions look so similar to 
the test particle ones, despite the action of such a diffusion 
process. A detailed analysis of the orbital structure in 3D 
models remains to be done, and in particular of the prop- 
erties of the vertical motion within regular and chaotic 
regions. 

11. Models versus observations 

Before concluding, we now present a selection of test par- 
ticle and TV-body velocity distributions yielding a good 
match to the observed one, confront the quasi-a;i(2) orbit 
and chaotic orbit interpretations of the Hercules stream, 
paying also attention to the case of the Hyades stream, 
and discuss how the models could be further improved. 

Beside the parameters in the initial conditions of the 
simulations, the free model parameters are Ro/Rolr, 
the velocity scale specified by Vo (defined as the local cir- 
cular velocity in the axisymmetric part <I>o of the poten- 
tial for the iV-body simulation), which should lie between 
180 and 230 kms^^ (e.g. Sackett |1997| ), and the veloc- 
ity of the Sun (wg, Ug) relative to the circular orbit in $0. 
A commonly adopted velocity reference in the Solar neigh- 
bourhood is the LSR, defined as the velocity of the most 
nearly circular closed orbit passing through the present lo- 
cation of the Sun according to Binney & Merrifield ( 19981 ). 
This definition, which is merely an attempt to generalise 
the circular LSR orbit of the axisymmetric case to non- 
axisymmetric potentials, is not always well adapted. The 
most reasonable LSR orbit candidates near the OLR of 
a barred potential indeed are the stable low-eccentricity 
a;i(l) and a:i(2) orbits, but some space positions near the 
OLR circle are visited by neither of these orbits in our 
models (see for example Ro/RohR = 1-0 and ip — 30° in 



Fig. |T^. However, for Ro/Rqlr ^ 1-0, there always exists 
a prominent peak of low-eccentricity quasi- 2:1(1) orbits in 
the model velocity distributions, which, as pointed out 
m Sect. |7[ not necessarily coincides with the trace of the 
non-resonant a;i(l) orbit when there is one. The maximum 
of this peak will therefore be taken as the model "LSR" 
and will be preferably associated to the Coma Berenices 
stream, which is the local maximum in the observed veloc- 
ity distribution that lies closest to the heliocentric velocity 
of the LSR ('1',m)lsr ~ (~5, 10) kms~^ derived from the 



Hipparcos data (Dehnen & Binney 1998). 

For Ro ^ -RoLR, the quasi- 2:1(1) peak is always close to 
the circular orbit of the axisymmetrised potential, except 
near the OLR radius and (p — 90° where it reaches a 
maximum positive w-offset of ^ 0.05uo for all explored bar 
strengths. Under the above circumstances and for realistic 
space positions, the azimuthal velocity of the Sun should 
exceed the circular velocity by 5 — 10 kms~^. 

The selected model velocity distributions are displayed 
in Fig. The distributions are derived according to the 
same procedures as described in Sects. ^ and Frame (a) 
shows a case where the Hercules-like stream is induced ex- 
clusively by chaotic orbits and peaks inside the H12 con- 
tour, illustrating the fact that chaotic modes not neces- 
sarily only occur in the hot orbit region. Here the Hyades 
stream coincides with a chaotic overdensity associated 
with a narrow and low-H chaotic breach roughly along the 
OLR curve, i.e. an interpretation similar to the one pro- 
posed in Sect. ^ for the u < extension of the LSR mode. 
Frame (c) gives a case where the Hercules-like stream now 
falls entirely in the hot orbit region and where the Hyades 
stream has the same chaotic origin as in frame (a). 

Frame (e), derived from the A^-body simulation and 
presenting a larger velocity dispersion, provides a remark- 
able example of a Hercules-like stream sustained exclu- 
sively by quasi-a;i(2) orbits. The test particle simula- 
tions develop quasi- 2:1(2) modes which cannot be as easily 
matched to the Hercules stream in our scaling procedure, 
generally falling right between this stream and the Hyades 
stream. This can be explained by the different local slope 
of the circular velocity Vc in the A^-body and the test par- 
ticle models. As explained by D2000 in terms of orbital 
frequencies, the separation between the quasi- 2:1(1) and 
the quasi-2i(2) modes increases with dwc/di?(i?o)- Since 
the average Af-body model has a slightly raising rotation 
curve near the OLR radius (Fig. |l6|d), its circular velocity 
gradient is larger than for the flat rotation curve test par- 
ticle simulations and thus the quasi- 2:1(2) mode is found at 
higher asymmetric drift relative to the quasi- 2:1(1) mode. 
However, the fact that observations support a rather gen- 
tly declining rotation curve at Ro and that a large inclina- 
tion angle of the bar is needed {(p ^ 40° for Ro > Rqlr) 
are arguments against the quasi- 2:1 (2) interpretation of 
the Hercules stream. On the other hand, the displace- 
ment of the quasi-xi(2) peak towards the H12 contour 
with increasing integration time mentioned in Sect. for 
F = 0.10 and Ro ~ i?oLR, may reinforce this interpreta- 
tion. 
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F = 0.10; F = 0.10; F = 0.15; F = 0.1.5; Af-body; 

JJo/iJoLR = 1.0.5; = 20°; iJo /iJoLR = 1.00; = 30° ; JJo /iJoLR = 1-05; = 10° ; iJo/iJoLR = 1.05; = 40°; iJo /iJoLR = 0.9; = 50° ; 

ti„ = 200; tia ^7; Us ^ -5. Vo = 210; t>, = 9; tl, = -15. t>„ = 180; ti^ = 5; tl^ = -7. Vo = 180; = 9; = -7. Vo = 200; t>, = 22; tl^ = 0. 
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Fig. 20. Selection of scaled velocity distributions from the test particle (frame a to d) and A'^-body (frame e) simulations, with 
the various parameters indicated on the top of the frames and the velocity origin at the adopted Solar motion. The velocity 
window and the velocity contours are the same as in Fig. |l[ The filled circles represent the mean velocities of the streams 
listed in Table |l|, excluding the Arcturus stream. All distributions from the test particle simulations are time averages over 
55 < t/tb < 65. The //-contours and the resonance curves are as in Fig. |^ and are not plotted for the A'^-body model because 
of the time delay problem discussed in Sect, hol 



Frame (d) is an example with two distinct low angular 
momentum peaks, the most negative v one being related 
to chaotic hot orbits and the other one to quasi-xi(2) or- 
bits. It would be interesting to check whether a sufhciently 
negative dvc/dR{Ro) is able to shift the quasi-xi(2) mode 
more towards the true location of the Hyades stream and 
thus yield a model velocity distribution with a better over- 
all match to the observed one. Note that the chaotic orbit 
mode will not necessarily be shifted as the quasi-a;i(2) 
mode, because its location in the u — v plane does not 
actually depend on the local slope of the circular velocity, 
but rather on the difference of $cff between the current 
space position and the Lagrangian point ii/2, which de- 
termines the u — V location of the i/12 contouiQ Finally, 
frame (b) displays a surprising case where the velocity dis- 
tribution in the quasi-a;i(2) region of the u — v plane (see 
Fig. H) seems to have split into two peaks coinciding well 
with the locations of the Hercules and Hyades streams, 
i.e. both these streams have a quasi-a;i(2) origin. However, 
this is likely to be a transitory situation resulting from the 
unachieved phase mixing near Ro = Rq^r (see Sect. 

These examples illustrate the variety of possible inter- 
pretations for the Hercules and Hyades streams, and it is 
very hard at this stage to decide with certainty which 
one is the most appropriate. The splitting of the LSR 
mode into the Pleiades, Coma Berenices, Sirius and other 
streams is probably not related to the bar itself and has 
a more local origin, like for instance the effect of time 
dependent spiral arms. 

In particular, at Rolr ~ 1-1 and 95 — 25°, where the low 
angular momentum mode has a chaotic origin, the v squash- 
ing of the bimodality reported by D2000 when decreasing his 
rotation curve parameter /3 is perhaps not the consequence of 
a local change of the circular velocity slope, but of an implied 
lower relative value of 'I'cff (/'i/2)- 



12. Conclusion 

The Galactic bar induces a characteristic splitting of the 
disc phase space into regular and chaotic orbit regions, 
with the latter regions owing only to the non-axisymmetric 
part of the potential in the limit of no vertical motion. In 
this paper, we have isolated these two kind of regions, 
as well as the quasi-periodic orbit sub-regions inside the 
regular regions associated with the stable a;i(l) and xi{2) 
periodic orbits respectively, within the same analytical 2D 
rotating barred potential with flat azimuthally averaged 
rotation curve as in D2000. We then have run test parti- 
cle simulations in this potential and a more realistic self- 
consistent 3D A^-body simulation to find out how the disc 
distribution function outside the bar region relates to such 
a phase space splitting and in particular how chaos may 
explain features in the Solar neighbourhood stellar kine- 
matics like the Hercules stream. 

Beside the bar strength, the regular versus chaotic 
splitting of phase space, investigated via the largest 
Liapunov exponent, is mainly determined by the value 
of the Hamiltonian H (or Jacobi's integral) and by the 
bar related resonances. In two dimensions and at fixed 
space position, the constant-/? contours in the galacto- 
centric u — v velocity plane are circles centred on {v,u) = 
{Roflp,0) and of radius \/'2{H — ^°ff), where Ro is the 
galactocentric distance, Vlp the rotation frequency of the 
bar and the local effective potential. The fraction of 
chaotic orbits increases with H and there is a sharp av- 
erage transition from regular to chaotic behaviour in the 
u — V plane when crossing the contour corresponding to 
the effective potential at the saddle Lagrangian points, 
H12 = $off(ii/2)j which generally intersects the w-axis at 
lower velocity than the circular orbit in the axisymmetric 
part of the potential. At H < H12, the orbits are rather 
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regular, while at if > H12, which defines the category of 
hot orbits susceptible to cross the corotation radius, they 
are rather chaotic. 

The resonances, on the other hand, generate an al- 
ternation of regular and chaotic orbit arcs in the velocity 
plane which, contrary to the low-t; part of the _ff-contours, 
are opened towards lower angular momentum. At bar in- 
clination angles (/? = and Lp — 90°, the maxima or min- 
ima of these stochasticity waves are in phase with the 
resonance curves derived from the axisymmetric limit and 
the arcs are symmetric in u, reflecting the four- fold sym- 
metry of the potential. At intermediate angles, these ex- 
trema become offset with respect to the resonance curves 
and the u-symmetry breaks. In particular, at i? ^ i?oLR 
and ^ 30°, a prominent regular region of eccentric 
quasi- a;i(l) orbits extends well within the hot orbit re- 
gion at w ^ 0, while the u > counterpart of the OLR 
curve is surrounded by a wide chaotic region consistent 
with the location of the Hercules stream. 

For a moderate bar strength {F — 0.10), the low- 
eccentricity and non-resonant quasi-xi(2) orbit regions ex- 
ist only for -R/i?oLR ^1-1 and over an angle range around 
ip — 90° increasing towards smaller R. There is no such 
region near the default position considered in D2000, i.e. 
i?/i?OLR ~ 1-1 and if — 25°, compromising the quasi- a::i(2) 
orbit interpretation given by Dehnen for the Hercules- like 
mode occuring in his simulations at the most realistic po- 
sitions of the Sun relative to the bar. 

The test particle simulations, started from axisymmet- 
ric initial conditions and progressively exposed to the full 
rotating barred potential, reveal a decoupled evolution 
of the disc distribution function within the regular and 
chaotic phase space regions. In the regular regions, the 
phase space density after phase mixing is roughly the same 
as the initial one, whereas in the chaotic regions, the par- 
ticles quickly evolve towards a uniform population of the 
easily available phase space volume via chaotic mixing, 
resulting in a substantial density re-adjustment. Because 
the space region within corotation represents a large initial 
reservoir of hot chaotic orbit particles which are spread 
throughout the disc by this process, yielding a net out- 
ward migration of such particles, the chaotic regions in 
the u — V plane outside corotation become more heavily 
crowded than the regular regions dX H > Hi2. In partic- 
ular, the wide and predominantly u > chaotic region 
mentioned above for realistic space positions of the Sun 
appears as an overdensity in the u — v distribution, provid- 
ing a coherent interpretation of the Hercules stream and 
explaining the u > property of this stream. According 
to this interpretation, the Hercules stream involves stars 
on essentially chaotic orbits which are forced to avoid the 
regular region at negative u. 

The time averaged disc u — v velocity distributions in- 
ferred from the A'^-body simulation are remarkably similar 
to those of the test particle simulations, despite the action 
of the transient spiral arms which allows at least some of 
the particles to diffuse from the regular to the chaotic re- 
gions and vice versa. At low eccentricity, the orbits are less 



sensitive to the inner features of the potential and the az- 
imuthal properties of the velocity distributions essentially 
align with the average local phase shift of the potential 
well relative to the bar major axis induced by the spiral 
arms. 

The velocity distributions may be very time dependent 
if for instance the bar has formed recently, because of the 
phase mixing occuring in the disc during at least ~ 10 bar 
rotations after the growth of the bar according to the test 
particle simulations. The u—v distributions in the A^-body 
simulation at fixed space position relative to the bar also 
display a strong temporal behaviour (see the mpeg movies 
at tittp : //www. mso . anu. edu. au/~fux/streains .html), as ex- 
pected from the presence of the transient spiral waves. 
However, since the simulation has been run for only about 
1.25 Gyr after the formation of the bar, phase mixing is 
still operating in the disc component, rendering difficult 
to disentangle from it the individual effects of such waves. 
The A^-body simulation also gives some insight on the con- 
sequences of evolving bar parameters: the slowly decreas- 
ing pattern speed of the bar mainly introduces a delayed 
response of the disc distribution function to the outward 
moving resonances, so that the velocity distributions at a 
given time rather reflect a higher value of R/Rqlr than 
the true instantaneous one when compared with the con- 
stant Qp test particle simulations. For completeness, one 
should mention that other perturbations than the bar and 
the spiral arms may provide alternative explanations of 
the local stellar streams, like for example the interactions 
of the Milky Way with its satellite companions. 

Finally, the process of chaotic mixing, combined with 
the possible stellar exchanges between the regular and 
chaotic phase space regions resulting from the diffusion 
of stars by transient spiral arms or molecular clouds, may 
provide an new and efficient way of heating galactic discs 
which remains to be explored. 
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